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Resume 

The effectiveness of projection methods for solving systems of linear inequalities is investi- 
gated. It is shown that they have a computational advantage over some alternatives and that 
this makes them successful in real-world applications. This is supported by experimental evi- 
dence provided in this paper on problems of various sizes (up to tens of thousands of unknowns 
satisfying up to hundreds of thousands of constraints) and by a discussion of the demonstrated 
efficacy of projection methods in numerous scientific publications and commercial patents (dea- 
ling with problems that can have over a billion unknowns and a similar number of constraints). 

1 Introduction 

Projection methods were first used to solve systems of linear equations in Euclidean spaces in 
the 1930s |3TJ|55] and were subsequently extended to systems of linear inequalities in [T|, l62l l63] . 
The basic step in these early algorithms consists of a projection onto an affine subspace or a 
half-space. Modern projection methods arc much more sophisticated [71 [51 [TUl [T71 [Ml [551 
[Ml 1321 HSl [57] and they can solve the general convex feasibility problem of finding a point in 
the intersection of a family of closed convex sets in a Hilbert space. In such formulations, each 
set can be specified in various forms, e.g., as the fixed point set of a noncxpansive operator, the 
set of zeros of a maximal monotone operator, the set of solutions to a convex inequality, or the 
set of solutions to an equilibrium problem. Projection methods can have various algorithmic 
structures (some of which are particularly suitable for parallel computing) and they also possess 
desirable convergence properties and good initial behavior patterns [8l[26l[33l|34l|35l[5Tl[67j. The 
main advantage of projection methods, which makes them successful in real-world applications, 
is computational. They commonly have the ability to handle huge-size problems of dimensions 
beyond which more sophisticated methods cease to be efficient or even applicable due to memory 
requirements. This is so because the building bricks of a projection algorithm are the projections 
onto the given individual sets, which are assumed to be easy to perform, and because the 
algorithmic structure is either sequential or simultaneous, or in-between, as in the block-iterative 
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projection metliods or in the more recently invented string-averaging projection methods. The 
number of sets used simultaneously in each iteration in block-iterative methods and the number 
and lengths of strings used in each iteration in string-averaging methods are variable, which 
provides great flexibility in matching the implementation of the algorithm with the parallel 
architecture at hand; for block-iterative methods see, e.g., Pl[TU l [TB l [ ^ [55 1 [iD l li l H5 1 [57 1 ISil 
and for string-averaging methods see, e.g., [T51 [5^ 1 [S5] . 

The convex feasibility formalism is at the core of the modeling of many problems in various 
areas of mathematics and the physical sciences ; see [32l[33] and references therein. Over the past 
four decades, it has been used to model significant real- world problems in sensor networks |14j . 
in radiation therapy treatment planning |2H I52j , in resolution enhancement [27j , in wavelet- 
based denoising [3D], in antenna design [55], in computerized tomography [5T], in materials 
science |56| , in watermarking [58j , in data compression [60| , in demosaicking |61| , in magnetic 
resonance imaging |69j , in holography |70j , in color imaging |71j , in optics and neural networks 
[72], in graph matching |73j and in adaptive filtering [75], to name but a few. In these - and 
numerous other - problems, projection methods have been used to solve the underlying convex 
feasibility problems. 

We focus on the important subclass of convex feasibility problems in which finitely many sets 
are given and each of them is specified by a linear equality or inequality in the Euclidean space 
M^. For such problems, which arise in many important applications [32l [511 [52] , alternatives to 
projection methods are available (see, e.g., [3l[48] and the references therein), and it is therefore 
legitimate to ask whether projection methods are competitive. 

In this paper we address this question and show that projection methods are indeed very 
competitive in the environment of linear inequality constraints. In Section [2] we discuss their 
comparative performance for four different kinds of problems. In Section [3] we give some 
examples of their use in real-world applications from the research and the patent literature. 
Finally, we present our conclusions. 

2 Comparisons 

2.1 Examples of 2-set feasibility problems 

In a recent paper [48], the author asks in the title : "How good are projection methods for 
convex feasibility problems ?" and immediately (in the Abstract) states that : 

"Unfortunately, particularly given the large literature which might make one think 
otherwise, numerical tests indicate that in general none of the variants [of projection 
methods for solving convex feasibility problems] considered are especially effective 
or competitive with more sophisticated alternatives." 

As indicated in the Introduction, projection methods have been used to solve highly nonli- 
near complex problems involving a very large number of sets. Therefore, results based on the 
geometrically simple 2-set problems of |48j are vastly insufficient to draw general conclusions. 
In addition, we show in this subsection that the experiments reported in [48j use suboptimal 
versions of projection methods, which further questions the justification of the above-quoted 
general conclusion as to their effectiveness. 

The numerical experiments provided in |48j focus exclusively on the problem of solving a 
linear system of equations under a box constraint, namely 



find X eR- 



N 



such that 



Ax = 6, 

N 
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where A G M.^^^'^ [M < N) has full rank, b E R^^, and the problem is assumed to be feasible. 
We show that, even in this basic setting, projection algorithms implemented with standard 
relaxation strategies perform much better than indicated by the results in |48| . 

Let us denote by Pi and P2 the projection operators onto the closed affine subspace Si = 
{x e M.^ I Ax ~ 6} and the closed convex set S2 = Xj^j^ [ci, di], respectively. The first operator 
is defined by 

Pi:x^x-A'^{AA'^)-\Ax-b), (2) 

where A'^ denotes the transpose of A. This transformation can be implemented in various 
fashions. For instance, in many signal and image processing problems, the matrix A is block- 
circulant and hence diagonalized by the discrete Fourier transform operator, which leads to a 
very efficient implementation of Pi [J. Here, we adopt a QR decomposition approach. Let 



Pii 




(3) 



be the QR decomposition of A'^ , where Rn is an M x AI invcrtiblc upper triangular matrix 
Then (HI) yields 

Pi: x^x-Qn{Rny\Ax-b). (4) 

On the other hand, the projection P2X ~ {T^i)i<i<N of a vector x ~ {xi)i<i<N onto ^2 is 
obtained through a simple clipping of its components, i.e., for every i £ {!,..., TV}, tt; ~ 
min{max{a;,;, Ci}, di}. 

Two standard projection methods to solve ([1]) are the alternating projection method 

x^o^eM^ and (Vn e N) = x^"' + A„(PiP2x(") - x^")) (5) 

and the parallel projection method 

G and (Vn G N) x'-"+'^ = + A„ ( ^ x^") ) , (6) 



where (A„)„gN is a sequence of strictly positive relaxation parameters. If A„ = 1 in ([5|), we 
obtain the popular Projection Onto Convex Sets (FOGS) algorithm [321 174] : 

a;(°)GR^ and (Vn G N) = PiPax^"' . (7) 

The convergence of any sequence (x^"^)„gN thus constructed to a point in S'inS'2 was established 
in [15]. On the other hand, if A„ = 1 in ([6]) , we obtain the Parallel Projection Method (PFM) : 

P, -r(") -I- P„T-(") 

xt") G and (Vn G N) x<-^+'^ = ^ ^ ' . (8) 

The convergence of any sequence (x*^"^)„gN thus constructed to a point in Si n S2 was es- 
tablished in [S], see also [5]. In [IS], (O and ^ are used, together with variants featuring a 
construction of A„ at iteration n resulting from a line search procedure and without closed-form 
expression. However, as the numerical results of [15] show, these relaxation schemes do not lead 
to significantly better convergence profiles than those obtained with the unrelaxed algorithms 
POCS ([7]) and PFM ([8]). In addition, nothing is said regarding the convergence of ([5]) and ([6]) 
with such relaxation schemes. 

The potentially slow convergence of projections methods has long been recognized [38 ] [50 ] [62] 
and remedies have been proposed to address this problem in the form of adapted relaxation 
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Fig. 1: Average performance of the algorithms when M x N = 600 x 1000. 



strategies tliat guarantee convergence. In tlie case of ([5]), it was shown in [TU] that any sequence 
generated by the Extrapolated Alternating Projection Method (EAPM) 



c(°) e Si and (Vn £ N) 



.("+1) 



r(«) ||2 



where < p <2 and A„ = { ||FiP2a;(") - 

.1, 



if ^ ^2, 
if a;(") e 5*2, 



(9) 



produces a fast algorithm that converges to a solution to ([T]). This type of extrapolation scheme, 
which exploits the fact that 5*1 is an affine subspace, actually goes back to the classical work 
of [SD] ■ It has been further investigated in [111 ISS] a-nd has been extended recently to a general 
block-iterative scheme in [TU]. Acceleration methods have also been devised for the parallel 
algorithm ([6]). Thus, the convergence of the sequence produced by the Extrapolated Parallel 
Projection Method (EPPM) 



and (Vn G 




Jn)\\2 



|P2a;('^ 



,(n)||2 



Pia;(") +P2x(") - 2x(")||2 



if ^ S"! n ^2, 
if a;(") e 5*1 n 52, 



(10) 



to a solution of ([T|) was established in [34 . This type of parallel extrapolated method goes 
back to [52 and [SS] , and it has been refined or generalized in several places (351 [HI IM]- 
In particular, it has been shown in numerical experiments to be much faster than unrelaxed 
projection algorithms in various types of problems ranging from numerical PDEs to image 
processing [35 l [46 l [66 l [67] . 

In Figure [H we compare the numerical performance of POCS ([7]), PPM ([8]), EAPM ([9]), 
and EPPM ^ for problems of size M x A = 600 x 1000. As in [TU] [SS] [35] , the performance 
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Fig. 2: Average performance of the algorithms when M xN = 600 x 1000 and condition 
numbers are around 3 x 10^. 



of the algorithms is measured by the decibel (dB) values of the normalized proximity function, 
which is evaluated at the nth iterate x*^"' by 

||Pix(o) -x(0)||2 + llPaxW -xWp j ■ ^ ' 

This comparison is relevant because the computational load of each iteration resides essentially 
in the computation of the projection onto 5*1 and it is therefore roughly the same for all four 
algorithms. The results are averaged over 20 runs of the algorithms initialized with x^^"^ = PiO 
and p = X = 1.9. In each run a matrix A G [—0.5, 0.5]*^^^ and a vector x € [0, 1]-^ are randomly 
generated. The vector b = Ax is then constructed so as to obtain a feasible problem using Ci = 
and = 1 in ll]). As in [48j and many other studies, we observe that POCS is faster than 
PPM. However, EPPM is faster than POCS and EAPM is clearly the best method : on the 
average, it is about 60 times faster than PPM, 30 times faster than POCS, and it achieves full 
convergence in just 7 iterations. In addition, convergence to a feasible solution is guaranteed 
by the theory and the expression of the extrapolation parameter Kn in Q is explicit and it 
requires no additional computation. It is argued in Section 5 of [48] that "there is a significant 
difference between random and real-life problems (similar observations have been made for 
linear equations, where random problems tend to be well-conditioned [Reference], and thus 
often easier to solve than those from applications)." Let us observe that random matrices do 
show up in many real-life problems, see |37U76| and the references therein. In addition, as shown 
in Figure [21 the qualitative behavior of the algorithms in the presence of poor conditioning is 
quite comparable to that observed in Figure [1] (for the experiments of Figure [21 the condition 
numbers vary from 3 x 10^ to 3.5 x 10''^). 

We have consistently observed this type of performance for problems of various sizes. For 
instance, we report in Figure [3] on the same experiment as above on problems of size M x N ~ 
3000 x 7000. Here EAPM is about 45 times faster than PPM, 22 times faster than POCS, and 
full convergence is achieved in just 5 iterations. 
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Fig. 3: Average performance of the algorithms when M x N = 3000 x 7000. 

These experiments indicate that the resuhs in [IS] on the speed of convergence of POCS 
and PPM (and the variants proposed there featuring modest speed-up factors and lacking a 
formal convergence analysis) correspond to a suboptimal implementation of projection methods 
and are not representative of their performance, since drastic improvements can be achieved by 
appropriate relaxations. 

2.2 An example from image representation 

The problem with the largest number of unknowns in the Nctlib/CUTEr LP problem set 
used in [48j has M x N ^ 6,330 x 22,275 and (according to the on-line attachment to [48]). 
for that problem, all methods discussed in |3H] need 42 seconds or more to reach the stopping 
tolerance on a 3.06 GHz Dell Precision 650 workstation. Wc found among the problems from 
applications that wc have been investigating one that is over an order of magnitude larger and 
for which the projection algorithm recommended in |28j required only 25 seconds on the average 
on an Intel Xcon 1.7 GHz processor, 1 Gbyte memory 32 bit workstation using the SNARK09 
programming system [41] . We now give a brief description of this problem. 

A J X J digitized image is one that is subdivided into square-shaped pixels within each of 
which the image value is uniform. Sometimes alternative representations of an image are super- 
ior. For example, in computerized tomography [51], we use the blob basis functions advocated 
by Lewitt [59j in some series expansion methods to reduce artifacts in the reconstruction. Such 
a reduction is due to the fact that blob basis functions are smoother than pixel basis functions. 

The contribution to the image value at the center of any of the M = pixels by any of 
the N blob basis functions is known from the geometry of the representations. If we are given a 
pixel image to start with and would like to find a good blob representation for it, the task is to 
find the weights x to be given to the blobs so that their combined contributions approximate 
the pixel values. In mathematical terms, this problem can be formulated as 

find X e such that c < Ax < d, (12) 

where the bounds c and d have to be tight to ensure a good approximation of the pixel image by 
the blob image. (The entries in the matrix A are the values of the various blobs at the centers 
of the various pixels.) 
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In the experiments reported in [28] M x N = 59,049 x 51,152. The algorithm that was found 
most efficacious among those tried is the projection method caUed CART3++ : the average (over 
40 instances of the problem) time required by CART3"'""'' to find a solution to (fT2|) was less 
than 25 seconds. 

The algorithm CART3++ belongs to a large family of projection methods that are usually 
referred to as algebraic reconstruction techniques (ART). These were first introduced to the 
tomographic image reconstruction literature in |47| ; for a recent discussion, see [511 Chapter 
11]. CART3"'"'', just like the closely related ART3+ that is used to solve the problem discussed 
in the next subsection, has an interesting mathematical property : provided that the set of 
feasible vectors satisfying the inequalities in p2p has a nonempty interior, both CART3++and 
ART3+ will find a feasible solution in a finite number of iterations [55]. This is achieved by 
appropriately controlling the sequence of relaxation parameters associated with the individual 
projections. 

2.3 An example from intensity-modulated radiation therapy planning 

The goal of intensity-modulated radiation therapy is to deliver sufficient doses to tumors 
to kill them, but without causing irreparable damage to critical organs. This requirement can 
be formulated as a linear feasibility problem of the kind shown in ()12|1 . The interpretation in 
this application is that each component of x is a to-be-determined strength of radiation to be 
delivered to the patient in separate beamlets, the components of Ax are the resulting doses 
at M points in the patient's body, and c and d are provided by the radiation oncologist as 
the desired limits on these doses. Two of the authors of the present paper (W. Chen and G.T. 
Herman) have been working in this area with D. Craft, T.M. Madden, K. Zhang and H.M. 
Kooy of the Department of Radiation Oncology, Massachusetts General Hospital and Harvard 
Medical School, and what follows in this subsection is an outcome of this collaboration. 

In the clinical case that we use as an example we have M x N = 302,491 x 13,734. The 
number of nonzero elements in A is 62,226,127, which is less than 1.5% of the total number 
of entries of A, an important consideration for the efficacy of projection methods for solving 
the problem. There is an additional technical consideration : since it is impossible to deliver 
negative radiation, each component of x has to be nonnegative, which results in an additional 
13,734 inequality constraints. As mentioned at the end of the last subsection, we use ART3+ 
[52] to solve this feasibility problem. 

In clinical applications, it is considered desirable to find multiple feasible points, each of 
which is optimal according to its own criterion. A typical optimization task is "find a feasible 
point that results in the smallest total dose delivered to the liver." The associated functional 
is a linear one : it is the sum of those components of Ax that are associated with points in the 
liver. Recognizing the speed by which ART3-I- finds a feasible point, we propose to apply it 
repeatedly, to solve the linear optimization problem 

find X G that minimizes x subject to c < Ax < d. (13) 

Our method solves this problem by turning the objective function into an additional constraint 
and solving 

find X e such that c < Ax < d and x < p (14) 

using ART3-I-. By reducing p using a bisection search until we obtain (within a prespecified 
tolerance) the lowest value possible for it, we get a good approximation to a solution of ([T3|). 
This whole process is called ART3+0. 

The task of minimizing a linear functional subject to linear inequality constraints is the 
well-known Linear Programming (LP) problem and several software packages are available for 
solving it, see, e.g., |3]. To compare the efficiency of our proposed procedure with currently 
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popular standard approaches, we applied them to the problem (|13p for a patient with pan- 
creatic cancer. We used all methods to find just a feasible point (No Task) and also for eight 
different LP tasks representing various linear optimization criteria. The three algorithms with 
which we compared ART3-I-0 were the self-dual interior point optimizer, the primal simplex 
optimizer and the dual simplex optimizer in the commercial software package MOSEK version 
5. The results are reported in Figure H) Typically, for each task, ART3-I-0 used about one to 
two minutes and the MOSEK algorithms needed one to several hours on an Intel Xeon 2.66 
GHz processor, 16 Gbyte memory, 64 bit workstation. It is also noteworthy that the memory 
requirements of the MOSEK algorithms were at least twelve times as large as that of ART3-t-0. 



2.4 Examples from computerized tomography 

Computerized tomography is the problem of recovering an image from its measured (and 
hence not strictly accurate) integrals along M lines [STj . If we assume that the recovered image 
will be represented as a linear combination of TV basis functions (see Subsection 12. 2p . then 
the task is to find the vector x the components of which are the weights to be given to the 
basis functions. Due to the linearity of integration and based on the knowledge of the basis 
functions, we can produce an M x matrix A such that Ax is approximately the vector b of 
measurements. Since it is not likely that there is an x such that Ax = b, it is reasonable to aim 
instead at finding an x that minimizes 

\\b - Ax\\^ + \\x\\\ (15) 

where cr G K indicates our confidence in our measurements. As explained in Section 11.3 of |51| . 
this sought-after x is in fact the x part of the minimum norm solution of the consistent system 
of equations 

[U a A] ab, (16) 

ixi 



8 



where U is the M x M identity matrix. In the same section there is a derivation of a variant of 
ART that converges to the sought-after x, given by : 

it^"-* is the A/-dimensional zero vector, 
is the A^-dimcnsional zero vector, 

^(n+l) ^ ^in) ^^^^^.^^ 

with 



a I 



(b, -a7xH)_„(") 
7n=A ^"^^"i |. ^" , (18) 



where, for n S N, j„ = (n mod M) + 1, for 1 < j < M, Cj is the M-dimensional vector whose 
jth component is 1 and whose other components are 0, aj is the jth row of A and bj is the 
jth component of 6, and < A < 2. Recognizing that in one iterative step only one row of the 
matrix is needed and that in computerized tomography most entries of each row are zero, we 
see that an iterative step can be carried out very rapidly, provided that we have access to the 
locations and the values of the nonzero entries. If the memory of the computer is large enough, 
this can be accommodated by storing ^ in a row-by-row sparse representation, otherwise the 
locations and values of the nonzero entries can be generated within each iterative step by some 
rapid mechanism, such as the digital difference analyzer explained, e.g., in Section 4.6 of |51| . 

In Section 5.8 of [51j there is an exact specification of the so-called standard projection data 
that are used to evaluate various reconstruction algorithms in that book, the number of lines 
used in the standard projection data is M = 223,744. In the evaluations based on the standard 
projection data that are reported in [51] for reconstruction algorithms that use blob basis 
functions, the number of blobs used is TV = 51,152. The first experiment on which we report 
in this subsection used exactly the same arrangement. (For the experiments in this subsection, 
the input data were created and outputs were analyzed and illustrated using SNARK09 [¥!].) 

In this experiment we applied the ART algorithm of (flT)) and (fT8| with ct = 5 and A = 0.05 
to the standard projection data. In Figure Ela) we show the behavior of the objective function 
(jl5p as a function of iteration cycles (an iteration cycle is defined to be M iterations). It can 
be observed that the initial decrease in the objective function is very rapid. 

This desirable initial behavior is even more noticeable when we evaluate the algorithm not 
from the purely mathematical point of view of how well the objective function is reduced, but 
rather from the application point of view of how good are the reconstructed images. For this 
purpose, we report on the normalized mean absolute picture distance measure, as defined in 
|51j . To define this measure we need a J x J digitization of the test phantom for which the data 
used in the reconstruction were collected ; such a digitization for the phantom we used is shown 
in Figure [61(a) . In our definition of the measure we use and sl"2 to denote the densities of 
the vth pixel of the uth row of the digitized test phantom and of the reconstruction (which is 
obtained from the vector x^"-* of blob coefficients) , respectively. We define the distance measure 
as 

J J 



^(n) _ u=l v=l 



J J 



(19) 



E E 1^"^^ 

u—l v—1 



In Figure [^b) we plot r*^"^ for this experiment. It is seen that its minimum is reached at 
the seventh iteration cycle, i.e., when n = 7M. This reflects the fact that the minimization 
objective ()15p does not (and, in fact, it cannot in real applications where the phantom is not 
known to us) fully describe the application objective. For this reason it is standard practice 
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(a) Plot of the objective function (|15p . 
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(b) Plot of the distance measure (|19|) . 
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Fig. 5: Image reconstruction by ART when M x N = 223,744 x 51,152. 
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Fig. 6: Displays of a 243 x 243 digitized phantom (a) and of an ART reconstruction 
when M xN = 223,744 x 51,152 (b). 

in tomography |51| to stop the iterative process after a few iteration cycles and use the result 
at that time as the reconstruction. The digitization obtained from x'-'^^^'^ produced by this 
experiment is shown in Figure [6l^b) . The reconstruction is not perfect (as indeed it cannot 
possibly be since the measured data are only approximations of the line integrals assumed by 
the mathematics), but important features of the phantom are identifiable in the reconstruction. 
This ART reconstruction was carried out in 38.4 seconds on an Intel Core 1.6 GHz processor, 
2 Gbytc memory, 32 bit laptop. 

We wanted to compare the time needed by ART with the time needed to solve the system 
(|16p of consistent equations for the same data by the current implementation of the interior point 
method of MOSEK version 5 [3j. Unfortunately this could not be done, because the memory 
requirements of the MOSEK software were too large for our laptop. So we attempted to use a 
much more powerful Intel Xeon 2.66 GHz processor, 16 Gbyte memory, 64 bit workstation, but 
even the 16 Gbyte memory was too small to handle this problem using the MOSEK software. 
The importance of this memory requirement issue for the subject matter of this paper cannot be 
overemphasized : problems that routinely arise in real applications can be handled by projection 
methods using inexpensive laptops, while "more sophisticated alternatives" fail to produce 
any results even on much more powerful workstations due to their much greater demands on 
computer memory. 

In order to be able to compare the efficiency of ART with that of the interior point method in 
MOSEK we had to reduce M and N to about a ninth of their previously-used sizes. Thus, in the 
second experiment on which wc now report AI x N = 24,880 x 5,711. For this smaller example 
we ran both ART and the interior point method in MOSEK (with its default parameters) on 
the Intel Xeon 2.66 GHz processor, 16 Gbyte memory, 64 bit workstation. In Figure [7| we plot 
both the objective function and the distance measure for both algorithms as a function of time. 
From the point of view of the objective function, MOSEK needed over 5000 seconds to reach a 
value as low as ART reached in 10 seconds. The advantage of ART is more pronounced when 
considering the picture distance measure : the optimal value is reached by ART at 1.7 seconds 
(when n = lAM) while the interior point method never reaches a distance value that is as low 
as that of ART and it needs approximately 5000 seconds to reach its lowest distance measure. 

Since both AI and N are about a ninth of their previous sizes, wc report in Figure [8] on 
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Fig. 7: Image reconstruction by ART and the interior point method in MOSEK when 
M X N = 24,880 x 5,711. 
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Fig. 8: Displays of an 81 x 81 digitized phantom (a) and of an ART reconstruction 
when M X N = 24,880 x 5,711 (b). 

the 81 X 81 digitizations of the phantom and of the reconstruction a;^^'^^^-'. These are clearly 
inferior to the images in Figure [51 demonstrating the medical necessity for the larger system of 
equations. 

3 Published and patented results 

3.1 Scientific publications 

Here we give a brief glimpse into some recently published results that show the efficacy 
of projection methods for some large problems. In the problems discussed in |40| . the number 
of unknowns was 59,049. In the examples given in [52] (a paper devoted to radiation therapy 
planning), problems of the form ()12p were considered with the number N of unknowns only 515 
but the number of pairs of constraints AI = 128,688. In four out of the six cases reported there, 
the projection method ART3+ [5^ found a feasible point in less than three seconds, and in the 
remaining two cases a feasible point was found in less than 34 seconds. These times are for a 
standard PC, using an Intel Xeon 1.7 GHz processor and 1 Gbyte memory. The problems in [40l 
[52] are small compared to some of the other applications for which projection methods have been 
successfully used. In [19j (a paper devoted to reconstruction from electron micrographs), there 
are examples in which 16,777,216 unknowns are to be recovered from 4,587,520 measurements 
(each giving an approximate linear equality) and others in which 884,436 unknowns are to be 
recovered from 92,160,000 measurements. Projection methods were used in [19] to handle such 
large problems in a reasonable time. 

In a recent paper |56] it is shown that a variant of ART can be used for crystal lattice 
orientation distribution function estimation from diffraction data. One of the problems discussed 
in [55] has 1,372,000,000 unknowns and the number of equations is potentially infinite. They 
are randomly generated and a projection step can be carried out as soon as a new equation 
is available (an ideal use of a sequential projection method of the row-action type, see j20j). 
The result reported in the paper for that problem is the one obtained after 1,000,000,000 such 
projection steps. 

As for all methodologies, projection methods arc not necessarily the approach of choice 
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in all applications. However, in important applications in biomedicine and image processing, 
projection methods work well and have been used successfully for a long time. For example, 
an important application of reconstruction from projections is electron microscopy and some 
of the leading groups in that field consider the projection method "ART with blobs" to be the 
method of choice, see [13]. A mathematical reason for this is that for such problems the angles 
between hyperplanes or half-spaces, represented by linear equalities or linear inequalities as in 
([1]) and (fT2l) , are in general large (in the sense that the cosine of the angle between the normals 
of two randomly chosen hyperplanes in the system to be solved is likely to be near zero) due 
to the high sparsity in each of the rows of the system matrix. 

3.2 Commercial patents 

There is hardly better evidence for the value of projection methods than the many pa- 
tents for commercial purposes that include them. Projection methods are used in com- 
mercial devices in many areas. Unfortunately, if a device is truly commercial, then the 
algorithm that is actually used in it is proprietary and usually not published. Many 
commercial emission tomography scanners use now some sort of iterative algorithms. A 
prime example is provided by the commercially-successful Philips Allegro scanners (see 
http ://www. healthcare. philips. com/main/products/ and [29j). In x-ray computerized tomo- 
graphy (CT), there are reports emanating from companies that sell such scanners indicating 
that variants of ART are used in heart imaging ; an example is presented in [54j . 

The first EMI (Electric & Musical Industries Ltd., London, England, UK) CT scanner, 
invented by G.N. Hounsfield [53], used a variant of ART. For this pioneering invention, Houns- 
field shared the Nobel Prize with A.M. Cormack in 1979. Thirty years later (on September 29, 
2009), a patent was issued to Philips (Koninklijke Philips Electronics N.V., Eindhoven, The 
Netherlands) for a "Method and device for the iterative reconstruction of cardiac images" [77] . 
The role of projection methods is demonstrated by the following quote from the "Summary of 
the Invention" included in the Patent Description : 

"The iterative reconstruction applied here may particularly be based on an Algebraic 
Reconstruction Technique (ART) (cf. R. Gordon, R. Bender, and G.T. Herman : "Al- 
gebraic reconstruction techniques (ART) for three-dimensional electron microscopy 
and x-ray photography", J. Theor. Biol., 29 :471-481, 1970) or on a Maximum Like- 
lihood (ML) algorithm (K. Lange and J. A. Fessler : "Globally convergent algorithms 
for maximum a posteriori transmission tomography" , IEEE Transactions on Image 
Processing, 4(10) :1430-1450, 1995), wherein each image update step uses the pro- 
jections of a selected subset, i.e., projections corresponding to a similar movement 
phase." 

4 Conclusion 

In this paper we have shown that, whether or not alternative methods are applicable, cor- 
rectly implemented projection methods arc very efficient for convex feasibility problems with 
linear inequality constraints, especially for those that are large, sparse, and originate from 
real-life applications. 
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